/*
Dieser Algorithmus berechnet das Ausgangsspektrum, von welchem dann das Maximum gesucht wird.
Eingangsvariablen/Zgr.:
        
        - Window_Re ... Realteil vom NDFT Window
        - Window_Im ... Imaginärteil vom NDFT Window
        - Re ... Realteil des Signals
        - Im ... Imaginaerteil des Signals
        - sizeValues ... Anzahl der Abtastpunkte (:=sizeof(Re) bzw. sizeof(Im) bzw. sizeof(tsampled) )
        - Nfreq ... Größe des Frequenzarrays fdetect (:= sizeof(freq) ) 
        - freq ... das Frequenzarrays fdetect 
        - tsampled ... der Zeitvektor
        - fDet ... Zgr auf den vorab reservierten Speicher für Ergebnisse
        - maxIter ... Durchläufe (:= wieviele f sollen detektiert werden) 
        - sinArray ... Vorberechneten Sinuswerte
        - cosArray ... Vorberechneten Kosinuswerte
        - P ... Zgr auf den vorab reservierten Speicher für das Spektrum  
        - P_clean ... Zgr auf das Filterspektrum
        - P_erg ... Zgr auf das Spektrum, welches bei der suche nach weiteren Frequenzen verwendet wird         

******************************************************************************************************
Status: noch keine Vereinfachungen, blanke Implementierung.
*******************************************************************************************************/


#include"tS.h"


void clean(double* Re,double* Im,int sizeValues,int NFreq,double* freq, double* tsampled,double* fDet,int maxIter, double * Window_Re, double * Window_Im,double* sinArray, double* cosArray, double* P_Re, double* P_Im, double* P_abs, double* Erg)
{
        double x,y;
        int start,ende,hilf;
        int NFreq_quart=(NFreq-1)/4;
        hilf=NFreq_quart*2+1;

        start=NFreq_quart;
        ende=start+hilf;
        NDFT(Re,Im,sizeValues,freq,tsampled,start,ende,P_Re,P_Im,P_abs,Erg,sinArray,cosArray);


        fDet[0]=Erg[1];
        fDet[0]++;
        double Ampl_Re=P_Re[(int)Erg[1]];
        double Ampl_Im=P_Im[(int)Erg[1]];


        for(int b=0;b<(maxIter-1);b++)
        {
                start=2*NFreq_quart-Erg[1];  //Update der Grenzen
                ende=start+hilf;
                int index=0;

                omp_set_num_threads(CORES);
                #pragma omp parallel private(index,x,y)
                {

                        #pragma omp for schedule(static)
                        for(int o=start; o<=ende; o++)
                        {
                                x= Ampl_Re * Window_Re[o] - Ampl_Im * Window_Im[o] ;
                                y= Window_Im[o] * Ampl_Re + Ampl_Im * Window_Re[o] ;
                                P_Re[o-start]-=x;
                                P_Im[o-start]-=y;                       //Subtraktion;kein gain-Faktor                  
                        }

                        #pragma omp barrier

                        #pragma omp for schedule(static)
                        for(int o=0;o<hilf;o++)                 //bereinigtes Spektrum:abs()
                        {
                                P_abs[o]=sqrt(P_Re[o]*P_Re[o]+P_Im[o]*P_Im[o]);
                        }
                }
                peak(Erg,hilf,P_abs);
                Erg[1]++;
                fDet[b+1]=Erg[1];
        }


}

//EOF
